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Abstract. The theory of adiabatic invariants has a long history, and very important 
implications and applications in many different branches of physics, classically and 
quantally, but is rarely founded on rigorous results. Here we treat the general time- 
dependent one-dimensional harmonic oscillator, whose Newton equation q + ui 2 (t)q = 
cannot be solved in general. We follow the time-evolution of an initial ensemble of 
phase points with sharply defined energy Eq at time t = and calculate rigorously 
the distribution of energy E\ after time t = T, which is fully (all moments, 
including the variance /i 2 ) determined by the first moment E\. For example, [i 2 = 
E 2 [(Ei/E ) 2 — (uj(T)/uj(0)) 2 ]/2, and all higher even moments are powers of ^ 2 , whilst 
the odd ones vanish identically. This distribution function does not depend on any 
further details of the function oj(t) and is in this sense universal. In ideal adiabaticity 
Ex = u(T)Eq/l)(0), and the variance /i 2 is zero, whilst for finite T we calculate E\, 
and /j, 2 for the general case using exact WKB-theory to all orders. We prove that if 
uj{t) is of class C m (all derivatives up to and including the order m are continuous) 
/i oc y-( m + 1 ) j whilst for class C°° it is known to be exponential /i oc exp(— aT). 
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1. Introduction 



In time-independent (autonomous) Hamiltonian systems the total energy of the system 
is conserved by construction, i.e. due to the Hamilton equations of motion, and the 
Liouville theorem applies, because the phase flow velocity vector field has vanishing 
divergence. In time-dependent (nonautonomous) Hamilton systems the total energy 
is not conserved, whilst the Liouville theorem still applies, the phase space volume is 
preserved under the phase flow. If (some parameter of) the Hamilton function varies 
in time, the energy of the system generally also changes. But, if the changing of the 
parameter is very slow, on the typical time scale T, there might be a quantity /, a 
function of the said parameter, of the energy E and of other dynamical quantities, 
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which is approximately conserved. It might be even exactly conserved if T — > oo, i.e. if 
the variation is infinitely slow, to which case we refer as the ideal adiabatic variation. 
Such a conserved quantity is called adiabatic invariant, and it plays an important role 
in the dynamical analysis of a long-time evolution of nonautonomous Hamilton systems. 

The theory of adiabatic invariants is aimed at finding the adiabatic invariants / and 
analyzing the error of its preservation at finite T. Namely, the statement of exactness 
of / is asymptotic in the sense that the conservation is exact in the limit T —>■ oo, whilst 
for finite T we see the deviation AI = If — 1$ of final value of // from its initial value 
Ii and would like to calculate A J. Thus for finite T the final values of / will have some 
distribution with nonvanishing variance. 

Here we just mention, as an example, that for one-dimensional harmonic oscillator 
it is known since Lorentz and Einstein pQ that the adiabatic invariant for T = oo is 
/ = E/uo, which is the ratio of the total energy E = E(t) and the frequency of the 
oscillator u(t), both being a function of time. Of course, 2ttI is exactly the area in the 
phase plane (q, p) enclosed by the energy contour of constant E. A general introductory 
account of the theory of adiabatic invariants can be found in 2j and references therein, 
especially [Sj, [I]. Another quite extensive excellent review is by Henrard [Sj, which 
includes rather complete list of relevant references to the original papers. 

However, in the literature this / and AI are not even precisely defined. It turns 
out that / must be generally considered as a function of the initial conditions, and then 
it turns out that it is conserved for some initial conditions but not for some others. 
As a consequence of that there is a considerable confusion about its meaning. Let us 
just mention the case of periodic parametric resonance (see section 5), with otherwise 
arbitrarily slowly changing u(t), in one-dimensional harmonic oscillator, where the total 
energy of the system can grow indefinitely for certain (almost all) initial conditions, and 
since u)(t) is bounded, I = E{t)/uj(t) simply cannot be conserved for the said initial 
conditions for all times, but only for sufficiently small times. In this work we consider 
/ as function of the initial conditions, and give a precise meaning to these and similar 
statements. 

Therefore to be on rigorous side we must carefully define what we mean by / and 
AI. This can be done by considering an ensemble of initial conditions at time t = just 
before the adiabatic process starts. Of course, there is a vast freedom in choosing such 
ensembles. Let us consider the one degree of freedom systems, which is the topic of this 
paper. If the initial conditions are on a closed contour /Q in the phase space at time 
t = 0, then at the end of the adiabatic process at time t = T they are also on a closed 
contour JCt which in general is different from /Co, but due to the Liouville theorem the 
area inside the contour is constant for all T. If /Co is a contour of constant energy at 
time t = 0, then JCt generally is not a contour of constant energy of the Hamiltonian 
at time T. The final energies of the system, depending on the initial conditions, are 
spread and thus distributed between some minimal and maximal energy, E min and E max , 
respectively. In the said case of periodic u>(t) and parameter resonance in a harmonic 
oscillator the contour fC T is squeezed in one direction and expanded in the other one 
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(the transformation is a linear map), and this contraction and expansion is exponential 
in time. Therefore, / = E(T)/u(T) can not be conserved; see section 5. 

Thus we must always study J as a function of initial conditions, by looking at the 
ensembles of initial conditions. There is a vast freedom in choosing such ensembles, but 
usually we do not know very much about the system except e.g. just the energy. 

Therefore in an integrable conservative Hamiltonian system the most natural and 
the most important choice is taking as the initial ensemble all phase points uniformly 
distributed on the initial iV-torus, i.e. uniform w.r.t. the angle variables. We call 
it uniform canonical ensemble of initial conditions. Such an ensemble has a sharply 
defined initial energy Eq. Then we let the system evolve in time, not necessarily slowly, 
and calculate the probability distribution P{E{) of the final energy Ei, or of other 
dynamical quantities. Typically E\ is distributed on an interval (E min , E max ), and 
P(Ei) is universal there, as it does not depend on any further properties of u(t) except 
for Ex. More general ensembles of initial energies w(E ) can be described in terms of 
the uniform canonical ensembles, as explained in section 6. 

To describe P(E 1 ) is in general a difficult problem, but in this work we confine 
ourselves to the one- dimensional general time- dependent harmonic oscillator, so N = 1, 
described by the Newton equation 

q + u 2 (t)q = (1) 

and work out rigorously P(Ei). Given the general dependence of the oscillator's 
frequency u(t) on time t the calculation of q(t) is already a very difficult, in fact 
unsolvable, problem. In the sense of mathematical physics is exactly equivalent to 
the one-dimensional stationary Schrodinger equation: the coordinate q appears instead 
of the probability amplitude ip, time t appears instead of the coordinate x and w 2 (t) plays 
the role of E — V(x) = energy - potential. If E is greater than any local maximum of 
V(x) then the scattering problem is equivalent to our ID harmonic oscillator problem. In 
this paper we solve the above stated problem for the general one- dimensional harmonic 
oscillator which is a follow up paper of our Letter 

In performing our analysis, we shall answer the questions as to when is / = 
E(T)/u(T) conserved, and if it is not conserved, what is the spread or variance /x 2 
of the energy, and the higher moments etc. Then (the not sharply defined) AI in 
the literature is proportional to fi. After performing the exact analysis, we provide a 
powerful technique based on the WKB method [7] to calculate /i 2 , and show that it gives 
exact leading asymptotic terms when T — > oo, and moreover, generally we can do the 
expansion to all orders, exactly. We treat several exactly solvable cases, and compare 
them with the WKB results, and finally prove the theorem as for how /i 2 behaves when 
u(t) is of class C m , which means having m continuous derivatives. 

We give a brief historical review of contributions to this field. After Einstein pQ, 
Kulsrud [8j was the first to show, using a WKB-type method, that for a finite T, I 
is preserved to all orders, for harmonic oscillator, if all derivatives of u vanish at the 
beginning and at the end of the time interval, whilst if there is a discontinuity in one of 
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the derivatives he estimated AI but did not give our explicit general expressions (37) and 
(38). Hertweck and Schliiter [§] did the same thing independently for a charged particle 
in slowly varying magnetic field for infinite time domain. Kruskal, as reported in jllj . 
and Lenard [TUj studied more general systems, whilst Gardner [TT] used the classical 
Hamiltonian perturbation theory. Courant and Snyder [TJ] have studied the stability of 
synchrotron and analyzed / employing the transition matrix. The interest then shifted to 
the infinite time domain. Littlewood [13J showed for the harmonic oscillator that if uj{t) 
is an analytic function, I is preserved to all orders of the adiabatic parameter e = 1/T. 
Kruskal ^3] developed the asymptotic theory of Hamiltonian and other systems with 
all solutions nearly periodic. Lewis [IS], using the Kruskal's method, discovered a 
connection between / of the 1-dim harmonic oscillator and another nonlinear differential 
equation. Later on Symon [IB] used the Lewis'es results to calculate the (canonical) 
ensemble average of the / and its variance, which is the analogue of our E\ and fi 2 . 
Finally, Knorr and Pfirsch ^7] proved AI oc exp(— const /e). Meyer [TH] . [T§] relaxed 
some conditions and calculated the constant const. Exponential preservation of / for 
an analytic u on (— oo, +oo) with constant limits at t — > ±oo, is thus well established 

m 



2. Transition map and general exact considerations 

We begin by defining the system by giving its Hamilton function H = H(q,p, t), whose 
numerical value E(t) at time t is precisely the total energy of the system at time t, and 
for the general time-dependent one- dimensional harmonic oscillator this is 

H = m + \ M ^ 2 ' (2 > 

where q,p,M,u are the coordinate, the momentum, the mass and the frequency of 
the linear oscillator, respectively. The dynamics is linear in q,p, as described by ([I}, 
but nonlinear as a function of oj(t) and therefore is subject to the nonlinear dynamical 
analysis. By using the index and 1 we denote the initial (t = to) and final (t = ti) 
value of the variables, and by T = ti — t we denote the length of the time interval of 
changing the parameters of the system. 

We consider the phase flow map (we shall call it transition map) 




Because equations of motion are linear in q and p, and since the system is Hamiltonian, 
$ is a linear area preserving map, that is, 



(4) 
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with det($) = ad — be = 1. Let E = H(q ,p ,t = t ) be the initial energy and 
Ei = H(q 1 ,p 1 ,t = ti) be the final energy, that is, 



E x = - 
2 



1 ( (cg + dp 



M 



+ Muf(aq + bp ) 



(5) 



Introducing the new coordinates, namely the action I = E/uj and the angle 



2E, 



MujI 



o / 

~2 cos 0, po = v 2ME sin < 



from (jSJ) we obtain 

-Ei = i?o( acos2 + /3sin 2 + 7 sin 20), 

where 



(6) 



(7) 



a 



+ a 2 % (3 = d 2 + ufM' 2 b 



2 ji -r 2 r 2 



M 2 a; 2 



7 = h abM—. 



Given the uniform probability distribution of initial angles equal to 1/(2-71"), which 
defines our initial uniform canonical ensemble at time t = t , we can now calculate the 
averages. Thus 



1 

2^ 



En 



En 



(a + /3). 



That yields E\ — E\ — E (5 cos 20 + 7 sin 20) and 



fi 2 = (Ex - E x 



— TP 2 
.2 _ ^0 fs2 



(S 2 + 7 2 ) 



where we have denoted 5 = (a — (3)/2. It follows from 
also in the form 

2 



/1 2 = (Ei - Ex 



El 



Ei 

Eq 



(9) 
(10) 

(jHJ) that we can write ^TQj) 

(11) 



As we shall see (see subsection 4.2), in an ideal adiabatic process [i = 0, and therefore 
E\ — Ex — luiE /ujo, and consequently P{E\) is a delta function, 



P(E 1 ) = 5(E 1 -u 1 E /u ). 



(12) 



Now we calculate higher moments of P{E\). Let us show that in general for 
arbitrary positive integer m 



and 



(Ei - Ei) 2m ~ x = 

(2m- 1)!! 



(Ei - Ei] 



2 m 



ml 



(Ei - Eif 



(13) 



(14) 



Energy evolution in time- dependent harmonic oscillator 



6 



It is easy to check the correctness of ()13|). so we prove only ()14|). Indeed, 



(Ex-E^ 2 ™ = £ 2m ( 7 sm20 + 5cos20) 2m = (15) 

2m 



k=0 

Note that for odd k = 21 - 1 



E ° m ( D (7 2m "^ fc sin 2m - fc 20 cos fc 20) . 



sin 2m - fe 20cos fe 20 = 0. (16) 
Therefore we can rewrite (|T3jl in the form 



2m\ 



(Ei - #i) 2m = El m J2 ( 2 ) (7 2m " 2 ^ 2 ' sin 2m " a 20 cos 2 * 20). 



(17) 



Using formula 2.512 of [20] we obtain 



1 /" ■ 2(m n nJ 2/o, „ (2/ - 1) ! ! (2m - 2/ - 1 ) ! ! 
— / sin 2(m " 20 cos 2 ' 20 ^ - v 



2vr y 2 m m! 
The substitution of this value into (fT7|) yields 



( 2 2 7) (^ (a - 2 ' - (») 
= ^|(^_g__) (20) 



2 m m! _ \l , 

- ^ - (lE^wy , (22) 

that is, the formula (|14|) holds. Thus 2m-th moment of P(E\) is equal to (2m — 
l)!!/i 2m /m!, and therefore, indeed, all moments of P(E\) are uniquely determined by 
the first moment E\. 

We now mention that the final energy distribution function written down as 



2tt . 



(23) 

<l>=<Pj(Ei) 



cannot be calculated analytically in a closed form in any useful way, because it boils 
down to finding the roots of a quartic polynomial, so we do not try to do that here, 
although numerically it shows interesting aspects (see figure 3). Obviously, P{E\) is in 
this sense universal, because it depends only on the average final energy E\ and the ratio 
uji/ujq of the final and initial frequencies, and does not depend otherwise on any details 
of uj(t). It has a finite interval as its support, between the lower limit E min and the upper 
limit E max , it is an even function there w.r.t. the mean value E\ = (E min + E max )/2, 
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and from E\ — E\ = E (5 cos 20 + 7 sin 20) it is easy to show that E min — E% — fiV2 
and E max = E\ + fiy/2, so that E max — E min = /ia/8 = 2E ^ 7 2 + 5 2 . At both extreme 
values it has an integrable singularity of the type 1 / \fx. This singularity stems from a 
projection of the final ensemble at t\ = T onto the curves of constant final energies E\ 
of H(q,p,ti). In between for every value of E\ = const = Ei(<f>), this equation has four 
solutions, namely 01,02,03,04, and thus we have to sum up all four contributions in 
the general formula ([23)1 . Of course, all that we say here for the distribution of energies 
Ei holds true also for the final action, the adiabatic invariant Ii = Ei/ui. It is perhaps 
worthwhile to mention that the moments of our distribution according to (jl4j) grow as 
2 m / ^iim, whilst e.g. in the Gaussian distribution they grow much faster, namely as 
2 m T(m + 1/2)/ \/7r, where T(x) denotes the gamma function. 

Expression (|TT|) is positive definite by definition and this leads to the first interesting 
conclusion: In full generality (no restrictions on the function u(t)\) we have always 
Ei > Equji/uoq and therefore the final value of the adiabatic invariant (for the average 
energy!) l x = E\/uj\ is always greater or equal to the initial value I = Eq/ujq. In 
other words, the value of the adiabatic invariant at the mean value of the energy never 
decreases, which is a kind of irreversibility statement. Moreover, it is conserved only for 
infinitely slow processes T = 00, which is an ideal adiabatic process, for which \x = 0. 
See section 4.3. For periodic processes u>i = u we see that always E\ > E , so the 
mean energy never decreases. See section 5. The other extreme opposite to T = 00 is 
the instantaneous (T = 0) jump where u switches to u)\ discontinuously, whilst q and 
p remain continuous, and this results in a = d = 1 and b = c = 0, and then we find 



Ei = ^(^ + 1), ^ 

Z UJ, 



2 



uj 2 



(24) 



y u L^O 

Below we shall treat the special case with u\ = 2lJq, and thus will find fi 2 /E 2 = 1/8 = 
0.125. 

Our general study now focuses on the calculation of the transition map (j4*)l. 
namely its matrix elements a, b, c, d. Starting from the Hamilton function (j2J) and its 
Newton equation (JTJ) we consider two linearly independent solutions ipiit) and ^if) and 
introduce the matrix 

m)-( Mt) Mt) ) (25) 

Consider a solution q(t) of (JTJ) such that 

q(t ) = go, q(to) = Po/M. (26) 
Because ipi and ip2 are linearly independent, we can look for q(t) in the form 

q(t)=AMt) + BMt)- (27) 
Then A and B are determined by 



* _1 (*b) q ° • (28) 
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Let q x = p 1 = Mq(t x ). Then from (H|-(|2BD we see that 




(29) 



We recognize the matrix on the right-hand side of (|2*9*j) as the transition map $, that is, 



3. Some exactly solvable special cases 

The theory of adiabatic invariants is rarely founded on rigorous results, therefore 
the study of exactly solvable cases is of fundamental importance, namely we can 
use them to test various analytic approximations and also the accuracy of numerical 
calculations. Because the ideal adiabatic processes are infinitely slow and refer to the 
limit T — > oo we must deal necessarily with the asymptotic behaviour of dynamical 
systems, which is difficult to approximate, because all simple minded perturbational 
and other approximation techniques fail (break down) after some finite time, and so 
typically give wrong predictions for the asymptotic behaviour in the limit T — > oo. The 
same difficulty occurs in numerical calculations. We shall see that the WKB methods 
[Zj can be successfully applied. In this section we deal with three different models for 
the function uj 2 (t). 

3.1. The linear model: class C° 

We assume that function uj 2 (i) is a piecewise linear function of the form 



Thus oj(t) has discontinuous first derivative at t — and t = T, and belongs to the class 
C°. Introducing the notation a = uj 2 , b = uj\ — uJq we obtain that on the interval (0, T) 
the equation has the form 




(30) 




(31) 



q+(a + -)q = 0. 

Two linear independent solutions of (|3*2*|) are given by the Airy functions: 



(32) 



bt + aT 

P/3TV3 



(33) 



and 



Bi{ 



bt + aT 

p/3^1/3 



(34) 
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The elements a, b, c, d of the matrix $^ are defined by the equation (J29|) and (J30j) . 
The exact analytic expression for 

E 



^ = {E l -E l y = ^(5 2 + 7 2 ) (35) 

is very complex, and we do not show it here. However, for — 1, u)\ — 2, Eq — 1, using 
the asymptotic expansion 10.4.60,62,64,66 of [21 j (pp. 448-449), we obtain the following 
approximation 

(Ei - E x y « ^ ^9 - 4 V2 cos(^^) J , (36) 
where we introduce the adiabatic parameter e, 

* = (37) 

We see in figure 1 that the exact expression (J33j) and its leading asymptotic 
approximation (|3*Bj) practically coincide, which demonstrates the power of the 
asymptotic expansion of the relevant expressions containing the Airy functions. Observe 
that the decay of fi 2 to zero as e — > is oscillatory but quadratic on the average, namely 
as y = jfge 2 , which means that \x goes to zero linearly with the adiabatic parameter 
e. This is always the case when uj{t) is of class C°. As we will see, in general, if 
u(t) is of class C m , then \i goes to zero oscillatory but in the mean as a power e m+1 . 
This theorem will be proven in subsection 4.2 using the exact formulation of the WKB 
method, applied to the relevant (but arbitrarily high) order. We shall see in subsection 
4.2, equation (jHEl) that the leading WKB term precisely reproduces the exact leading 
term in (J3BJ). In figure 2 we show the same thing as in figure 1 but on the larger scale 
of e. 

We now compute the final energy distribution function 



dEi 



1=1 

Note that 

d< f> = 1 (OO) 

dE x ew y } 

Using the substitution 

u = tan(0/2) (39) 
we obtain from ((7j) 

(1 + u 2 fE l = E (a(l - u 2 ) 2 + A(3u 2 + 4714(1 - u 2 )). (40) 



Resolving this equation with respect to u yields four branches of the function (f)(Ei). In 
figure 3 we show the emerging distribution function as an illustration. 
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Figure 1. (E\ — E\) 2 for < e < 0.05; the lines of the exact expression (|35() and 
the asymptotics (13611 practically coincide; the non-oscillating thin line is the parabola 
y i28 c • 




0.2 0.4 0.6 0.8 1 1.2 



Figure 2. [Ex — Ex) 2 for < e < 1.2; the lines of the exact expression l|35|) (full) and 
of the asymptotics (j36(l (dashed) practically coincide for e < 0.3; the non-oscillating 
thin line is the parabola y = yfgC 2 . One can show that /j 2 goes to 1/8 = 0.125 when 
e — > oo, which means T — > 0, which means the instantaneous jump of lu = 1 to 

LUx = y/2. 
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1. 4121. 4131. 4141. 4151. 4161. 417 

Figure 3. The distribution function P(E\) for the linear model of subsection 3.1, for 
E = l,ul = = 2,e = 0.01. 



3.2. The harmonic model: class C 1 

Now we consider the case when u(t) is of class C 1 , namely we use the specific model 



u\t) 



a if t < 

6 - (6 - 5) cos(^) < t < T 



(41) 



2b -a 



if t > T 



where a = u 2 , 2b — a = u\. Then the equation (JTJ) has the form 



Tit 



q+ (6 - (6 - 5) cos y—j jq = 

and the fundamental solutions are represented by the Mathieu functions: 

AbT 2 2{b-a)T 2 nt 



^i(t) = ce(- 



7T^ 



7T- 1 



2T' 



and 



,46T 2 2(6 -a)T 2 nt 



71- 



7[- 



2T' 



(42) 



(43) 



(44) 



In this model the first derivative of u(t) vanishes at t — and t = T and is continuous 
there, whilst the second derivative is discontinuous, therefore u{t) belongs to the class 
C\ 

Here we find that indeed \j? goes to zero oscillating but in the mean as e 4 , as shown 
in figure 4. The exact analytic asymptotic behaviour of the Mathieu functions needed 
for our calculations is not known, unlike for the Airy functions, but we compare the 
exact (numerical result) with the WKB approximation expounded in section 4, see the 
final general formula (J98|) of subsection 4.3. 
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Figure 4. We show, for the harmonic model of subsection 3.2, {E\ — Ei) 2 for 
< e < 0.1. The exact result is represented by the full line, whilst the dashed 
curve is the curve 0.056 e 4 (41 + 9 cos(^P)), obtained by the WKB method, equation 
(JSSJ of subsection 4.3. 



3. 3. The analytic model: class C c 
Consider now the case when 



- 2 W = Tj£- < 45 > 



9 + -r-^rv = o- (46) 



that is, ( [TJ is of the form 

1 + a e at 
1 + e c 

(Please do not confuse these a and a with the usage in the main text.) Then two 
fundamental solutions of (|4l)j) are represented using the hypergeometric function 2 Fi as 
follows: 



Mt) = e-» - i^, -1 + i^, 1 - -, -*«) (47) 

a a a a a 

and 

, / % ;+ „ / i i \fo> i i \/a 2 i „ , , . 

Mt) = e 2*i( — , - + — , 1 + -, -e at ). (48) 

a a a a a 

Using (|37jl and (|4^|) we compute the matrix $ defined by (JHUj) with to = — ^ and 
t\ = — , so T = ti — £ = 77 (note that the most essential change of the function u(t) 
is on the interval (—-,-))• For the case a = 2 the graph of (E\ — Ei) 2 as the function 
of e = 1/T is displayed in figure 5 as continuous line. The approximation to the graph, 
y = 4.174e~ 634//e , computed using the least square method, is displayed in the figure as 
dashed line. We clearly see that the behaviour of \j? is exponential. 
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0.035 0.04 0.045 0.05 



Figure 5. The variance (E\ ~ Ei) 2 of the energy for the analytic model (|45|l . for 
0.03 < e < 0.05. The dashed curve is approximation y = 4.174e _0 634 / e . 



4. Application of the WKB method for the calculation of the transition 
map 

4-1. General and exact considerations 

In this section we proceed with the calculation of the transition map $ in the general 
case, and because (TJQ) is generally not solvable, we have ultimately to resort to some 
approximations. Since the adiabatic limit e — > is the asymptotic regime that we 
would like to understand, the application of the rigorous WKB theory (up to all orders) 
is most convenient, and usually it turns out that the leading asymptotic terms are well 
described by just the leading WKB terms if e is sufficiently small. In using the WKB 
method we refer to our work [Jj, where we have derived the explicit analytic expressions 
for all WKB orders in closed form, except for the exact rational coefficients, which can 
be easily obtained from a recurrence formula, i. 

We introduce re-scaled and dimensionless time A 

A = et, e = 1/T, (49) 

so that (Q) is transformed to the equation 

eV(A)+^ 2 (A)g(A) = 0. (50) 

By prime we denote the differentiation w.r.t. A. Let q+(\) and ?-(A) be two linearly 

| There is substantial literature on WKB method, which due to limited space cannot be reviewed here. 
But we should mention the classic works by N. Froman and P.O. Froman, who have found a number 
of interesting relationships, e.g. a relation between the even and odd order terms |22|. although we do 
not use it here, so that our exposition is selfcontained. 
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vR = ( 1 (^) 

X y eMq' + (X) eMq'_(X) J 1 ' 

and taking into account that Ao = et , A x = eti we obtain for the matrix the 
expression 

*=[ a c dj =^a(A 1 )^ 1 (A ). (52) 

We now use the WKB method in order to obtain the coefficients a, b, c, d of the matrix 
$. To do so, we look for solution of ()5()j) in the form 

q(X) = wexp j-^a(A) j (53) 

where cx(A) is a complex function that satisfies the differential equation 

(a\X)) 2 + ea"(X) = -cu 2 (X) (54) 
and w is some constant with dimension of length. The WKB expansion for the phase is 

oo 

<7(A) = £Va fc (A). (55) 

fc=0 

Substituting (J55|) into (J53j) and comparing like powers of e gives the recursion relation 

n-l 



^ 2 (A), ^=-^(i>^- fc +<-i)- ( 5e ) 



/2 

fc=i 



Here we apply our WKB notation and formalism from our work [7j and we can choose 
°o+(A) = iw(A) or ctq_(A) = — iu;(A). That results in two linearly independent 
solutions of (J5Uj) given by the WKB expansions with the coefficients 

(T ,±(A) = ±i / u(x)dx, <ri,±(A) = -j-log ^ \ , (57) 



cr 2i± = ±- / —3 dx, ... (5£ 



i /" A 3c/(:r) 2 - 2u(x)u"(x) 

Since l<j(A) is a real function we deduce from (|5S|) that all functions cr' 2k+l are real and 
all functions o' 2k are pure imaginary and o~' 2k + = —a' 2k _, cr 2k+1 , = c^fc+i - where 
fc = 0,1,2,..., and thus we have cr+ = A(X) + LB (A), oi. = A(A) - LB (A) where 

^(A) = E^o e2fc+la Wi( A )' B ( X ) = -iE^o e2fca k+( A ) are both real quantities. 
Integration of the above equations yields 

a+ = r{X) +is(A), a_ = r(A) - is(A), (59) 
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where r(A) = J x A(x) dx, s(A) = f£ B(x) dx. Below we shall denote S\ = s(Ai). 
To simplify the expressions let us denote A = A(Aq), A\ = A(\\), B = B(\ ) and 
Br = B(Ai). 

Using this notation we find that the elements of the matrix $^ have the following 
form: 



e e 



.A,, mi) ( — j — £> cos^— 



e ^ . /si 
■sin — 



M£? \ < 

-DO 



(^Ai + BoSOsin^ 



si 



+ (A B 1 - A X B Q ) cos l^— 

Ai sin ^— j + i?i cos ^— 



e £ 



o 



Therefore for the quantities a and /3 defined by (JHJ) we obtain 



(60) 



a 



, ,2 E>2 



c^i ^4 sin ^) — -^o cos (^r)) + 



((A)Ai + S Si)sin (^) + 



(AqBi - AiB ) cos 



(61) 



e £ 



+ (X sin (-) + B L cos 



(62) 



We remind that det(\P) = 1. Computing now det(^) = ad — be using the expressions 
(HU) we find 

2 

(63) 



K*i) 
e « 



and therefore finally, exact to all orders, 



1 



a = — 



vBqBi 
1 



A sin ( — j — B cos ^— 



(64) 



Si 
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d 



+ (A B l - A 1 B ) cos 

A\ sin ^ — j + B\ cos — 



VBqBi 

Thus we obtain the final result for the expression Q, exact to all orders, 



a + (3 



BqBi 



sm 



2 (Si\ ( B\B\ 



+ UJ 1 J + cos 



^0-2 

a;. 



-o 



Sl \ ( 2 <4 A 2 A( 2A A 1 B B 1 



snr(-) 



+ 



COS 



2 / Sl N , A*B* 2A A 1 B B 



+ 



o 



+ 4f + 



+ 



sm | — 1 cos I — ) x 

,2 



2A B ^ + 2A 1 B 1 + A [AqA\ + B Bi) (4,£i - AM, 



UJn 



4-2. Leading asymptotic terms in the power expansion in terms of e 

So far the result is exact. Let us consider the first order WKB approximation, which is 
the generic case, that is 



(65) 



A(\)nea' 1)+ (\), B(\)^^±^-=u(\). 
Substituting these values of A(X) and B(X) into we find 



(66) 



a + p = 2— + 

UJ 



( 



£JlU/(V 2 

4ul 



cos 



2j^co(x)dx 



\ 



As we have shown above E\ = ^-^-Eq and 



2uqUJi 



(AEi) 2 _ l 



1} , "'(Ai) 2 



J 



(67) 



0(6 3 ). 



. Therefore 



El 
( 



Wiw'(Ao) 



V 



(68) 



( , )S , jA o y ) u/(A )u/(Ai) 



4^ 



+ 



^'(Ai) 2 



\ 



2, ,2 



+ 0(e 3 ). 



/ 



Substituting into the last formula u;(A) = y/1 + A we obtain exactly the approximation 
(pffij). Thus the WKB approach yields exactly the leading asymptotic term for E\ and 
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fi 2 for general u>(t), and not only that, but also all the higher power terms of e if desired. 
Below we shall explain how all the higher order WKB terms can be calculated using our 
closed form formula [Zj. Please note that the expression in the big brackets of (ffiH|) is 
positive definite (its minimal value is a complete square of a real quantity), as it must 
be. 

Now we consider the special cases where w'(Xq) and w'(Xi) vanish which implies 
that the leading order terms in the above two equations vanish. In such a special case 
we have to work out the higher order WKB terms. By this reasoning and calculation 
we prove that if u(t) is of class C m (all derivatives up to and including the order m are 
continuous) then [i oc e^ m+1 \ 

Therefore we now suppose that all derivatives at Ao and Ai vanish up to order 
(n-1) 

u/(A ) = . . . = w (n-1) (A ) = uA-^Ai) = 0, cu (n) (A )cj (n) (A 1 ) ^ 0. (69) 

Then 

*i(Ao) = ^(Ai) = • • • = <_!(Ao) = <_!(Ax) = 0, <(A ) ( r n (A 1 ) ^ 0. (70) 

There are two cases to be considered separately, namely, the case of odd and the case 
of even n. 

In the case n = 2k — 1 in order to find an approximation to fi 2 up to order e Ah ~ 2 
we expand A(X) and -B(A) as follows: 



2k 



A{\) = ^V^+CA) + 0(e 4fe+1 ), (71) 

u=k 

2fe-l 

B(A) = W (A)-i^ f % + (A) + 0(^. 



u=k 



We observe that the terms of the approximation (fTTj) of the order higher than e 2k are 
essential only in the following terms appearing in the brackets of the expression (|f)5|l : 

cos 2 m u 2 B 2 2 / Sl x _ sm 2 ( s -±)B 2 B 2 , s 

1 ' cos ( — ) B 2 , Ke> - , (72) 



\ e / uj 

sin - cos - I -2 1 9 ° + 2A 1 B 1 - 2 1 + 2 ? 1 
v e / \ e / \ u o 

Taking into account that 

2k-l 



B(X) 2 = u(X) 2 - 2w(A)i £ ^k + (A) + 0(e 4fc ), (73) 

u=k 

BlB\ = wgB? + ufB 2 - u 2 u 2 + 0(e 4fc ), (74) 
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and 

A(X)B(X) = uj(X)A(X) + 0(e 4k - 1 ) (75) 

we easily find: 

2A 1 B 1 -2^^ = 0(e 4k - 1 ), (77) 

yielding that the coefficient of sin (^) cos (7 1 ) m dZ2j) is of the order 0(e 4fc_1 ), and we 
find 



cos 2 (n) ulBl _/ Sl x _ sin 2 (sA BlB{ 

KeJ 1 + cos 2 ( — ) B\ + uy 2 1 = (7* 



c<7q v e / w, 

u\ sin 2 ( ^ + 



o 2 2 
20^ cos ^ — ^ 

2Jfc-l , ■■> s 

£ e 2 « ( -2o; 1 i< i+ (A 1 ) - 2i^< i+ (A ) ) + 0(e 4fc ~ 1 



Therefore the expression in the brackets of 1)65)1 is equal to 



2k-l , 2 \ 

2u l + e 2u -2o; 1 i4 U)+ (A 1 ) - 2i^< )+ (A ) + (79) 

Ak-2 ^l 2(T 2fc-l,+ ( A o) 2 , x 2 

e 2 1" °2fc-i,+( A oJ - 



2 ^l (T 2A : -l,+ (^o) Cr 2fc-l,+ (^l) . 2 / S l 



CJ 

Noting that 



(l_ 2s in 2 (^))) +0(e 4fc -^ 



1 _ ^(i + i 2 ^^ 2 ^^k±iM + ^±iM)) + 0(^) (80) 



x u=k 



we finally obtain 



(A^) 2 4fc _ 2 /<_i j+ (Ai) 2 , w 2 a 2fc _ li+ (A ) 2 



£g - £ I 2uJ + 2^ < 81 > 



W l < ^*-l,+ (*o) <7 2*-l,+(^l) / 2s l\\ /->, .I/.-J 



cos — + 0(e 



^0 V e 

The second case is when n is even, n = 2k, k positive integer, so that we assume 

now 
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2/v 



A(A) = £^V^ li+ (A) + /i. .f. I (82) 

2fe 

B(A) = w(A)-i£e^ + (A) + /i. .f. 
Then, similarly as above, using the equalities 

2Jfe-l 

5(A) 2 = ^(A) 2 - 2w(A)i £ e 2 <V^ + (A) - e* k (a' 2k>+ (\)) 2 + 0(e 4fe + 2 ), (83) 

M=fc 

BlB\ = a; 2 5 1 2 + c, 2 J B 2 - Wo 2 W 2 -4^i6 4 V^ i+ (Ao)^ fci+ (A 1 ) + 0( e 4fc+2 ), (84) 
A(A).B(A) = uj(X)A(X) + 0(e 4fc+1 ) (85) 

and 



we obtain 



.i, / ^ 2fc , + (Ao)) 2 (^, + ( A o))(^, + (Ai)) 

^0 ^0^1 




(A£i) 2 * k f *k, + M 2 , ^i 2 ^, + (Ao) 2 , S7 , 
~^""" C + (8T) 

^l<,+ (AoX,+ (Al) cQ Y2Si'| | + /4fc-l, 



We now observe that the expression 



E l V 2uJ l 2 ^o 

"'< +( y- +(Ai) cos(^)) + o(^'). 

coincides with (18 1|) when n = 2fc — 1 and with (|87|) when n = 2fc. If n-th derivative 
of cj(i) at Ao and Ai is nonzero, it means that u(t) is of class C m where m — n — 1, or 
n = m + 1. Then /i 2 is oscillating as e — > but in the mean goes to zero as a power 
/i 2 oc e 2 ( m+1 ) ; and therefore the general equation (J88j) proves our assertion. 
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4-3. Further simplifications of the general formula for the leading asymptotic term 

Now we want to simplify the expression (J88|) by expressing a' n in terms of u' n , using our 
explicit results in 0. 

Let M = U^N*, N is the set of non-negative integers. We define the map 
L : M -> N by 

L(v) = 1 ■ vx + 2 ■ v 2 + . . . + I ■ ui (89) 

and denote by L{y) = m the equation 

L{y) = 1 • Vx + 2 • u 2 + . . . + m ■ v m = m, (90) 

with m G N, v G M (there is one-to-one correspondence between the set {(u±, . . . , 

of solutions of (j9*Uj) and the set of partitions of m). For a vector ^ = . . . , i/j) G M we 

denote QM = (QTW)" 2 • • • (Q^)" 1 , IH = fi + • • • + *4 and let (i = 1, 1) 
be the vector {u\, . . . , z/j + 1, z/ i+1 — 1, . . . , i/j). It is shown in [7] that the functions a' m are 
of the form 



v:L(y)=m ^ 

where the coefficients Z/^ satisfy the recurrence relation 



U v = (92) 



1 TTTT , (4-LH-2|^|)^ 1 _ 1 ^„„^ ) ^ (^ + 1)^ } 
- 2^ ^ + 4 + 2 ' 

with f/g = — 1 and we put U a = if among the coordinates of the vector a there is a 
negative one. 

In our case Q(x) = —u 2 (x). Therefore in the case when 

u/(X) = ...= J n " l \\) = 0, U W(A) ^ (93) 
from (}9"Tj) we obtain 

< = ^kr, (94) 

Q 2 

where ZTj? is computed by (|9*2^1 and v — (0, . . . , 0, 1) is the solution of L{y) = n with all 
entries except of the last one equal to zero. From (|9*2"J) we obtain U& = 2^1^(1,0,. ..,0) = 
(| + |f/o) = -i^T yieding 

< = (95) 

If fl£3) holds then 

Q(») = -2uju {n) . (96) 
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Therefore, from (J95|) we obtain 

■n+l (n) 

Substituting this expression in (|K%|) we obtain the final general formula for the leading 
asymptotic term of /z 2 , namely 



(AE 



\2 



E 2 

J2n / .21. ,(n)\2 /, ,(")\2 , W (n) 



(9? 



o2n+l I 2(n+2) 1 \2n h ,2 n+3. ,n-l 




^0 

In the special case n = 1 we recover the formula (|68p. We emphasize that the above 
expression (j9~%|) is indeed evidently positive definite, which must always be the case for 
the variance. 

If the lowest nonvanishing derivative is the n-th derivative of u(t) at Ao and/or Ai, 
then it means that u{t) is of class C m where m — n — 1, or n — m + 1. Then /i 2 is 
oscillating as e — > but in the mean goes to zero as a power jj? oc e 2 ( m+1 ) ; and therefore 
the general equation proves our assertion. 

If o;(t) is an analytic function on the real time axis (— oo, +oo), the decay to zero is 
oscillating and on the average is exponential cx exp(— const /e) [T^j, [THJ , [IH]- This 
exponential smallness stems from the divergence of the relevant series and has been 
extensively studied in related works j2Hj, jSlj, [25] . [26| and [22], where the resummation 
techniques have been devised. 

Let us now summarize our results for the variance /i 2 as a function of uj{t) embodied 
in the exact general formulae and (jHSJ)- If ^{t) is analytic between to and ti > to 
then the equation (}9"%j) applies, and as we see /i 2 is dominated by the switch-on and 
switch-off events at to an d ti, respectively. However, the smaller the lowest nonvanishing 
derivatives of u(t) at the two points, the smaller will be the power law contribution. 
Indeed, if t and t x go to — oo and +oo, respectively, and u(t) is analytic on the entire 
interval, then the behaviour will become exponential at sufficiently large e > e c , whilst 
it will be oscillatory and in the mean a power law at small e < e c . For example, in figure 
5 we observe only exponential regime, because e c is so small. 

If there is any non-analyticity of u(t) on this interval, say at tj, then the calculation 
of /i 2 must be split into two intervals (to,tj) and (tj,ti), calculating the transition 
matrices for each interval separately, then looking at their product and calculating its 
matrix elements a, b, c, d, and then proceeding to calculate a + (3 and finally /x 2 , which 
will obey a power law like in (|98j) . 

In other words, if u(t) is analytic everywhere, and if the switch-on and switch- 
off events are completely eliminated, in the sense that they are infinitely smooth, the 
behaviour of /x 2 (e) is exponential, as explained above. In all other cases it is a power 
law. 
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5. Periodic uj(t) 

In case when u(t) is periodic with period r but otherwise completely general we can 
state some general rigorous results. Since the frequency Uq at time t and LO\ at time 
t\ = t + r are equal, we see from equation (fTT| due to positive defmiteness of /x 2 that E\ 
is always greater than Eq, that is, in a period r, or any integer multiple of it, T = nr, 
the mean energy E% never decreases. 

If we denote by $1 the transition map (J3J for one period of our periodic system, 
then the transition map $ n for an interval of exactly n periods of length r is simply a 
power of $1, namely 

$n = (99) 

If we use units such that Uq = uj\ = 1 and M — 1, then the mean energy from (|SJ) and 
(jHJ) can be expressed elegantly as the trace of the product of the transition matrix $ 
and its transpose $ T , namely 

E x = ^(a + (3) = |V + b 2 + c 2 + d 2 ) = ^Tr($$ T ). (100) 

If $ = $„ = $™, then we have to understand the behaviour of the above expression ()100j) 
as a function of the number of periods n. This is obviously dictated by the eigenvalues 
of the transition map $i, whose determinant is of course equal to 1, and its trace is 
denoted by S, which can be reduced to the diagonal form by a similarity transformation 

$i = WDW- 1 (101) 

where W is the transformation matrix and D is the diagonal matrix with the eigenvalues 
ei and e 2 = 1/ei, namely the solutions of the quadratic equation for e, 




-eS + l=0, 5 = Tr$x e x = l/e 2 = - ± \ - -1. (102) 



Therefore e±, e% — 1/ei are either real reciprocals (if \S\ > 2) or complex conjugated 
numbers on the unit circle (if \S\ < 2). The n-th power of $i can then be written as 

$ = = = WD n W~ 1 . (103) 

Therefore, the matrix elements a,b,c,d of $ n are bounded when IS*! < 2 and oscillate 
with n, while they increase exponentially when |>S| > 2. Indeed, if e± > 1, then the 
asymptotic behaviour of ()100j) is 

E 1 w KE e 2n , (104) 

where K is a constant determined by the matrix W in the transformation (jlOljl . and 
the variance at sufficiently large n goes asymptotically as 
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Therefore, in case when e\ > 1, we indeed find that the mean energy and with it 
the variance of the energy increase exponentially with time T = ra~, and since ui(t) 
is bounded, nothing is conserved. t\ is determined by S = Tr$i, and this in turn is 
determined by the specific properties of the system, that is by the nature of u(t) on the 
interval (0, r). 

So far we discussed the statistical properties of the energy distribution P(E n ) in 
a periodic system. Of course, the contour of the initial uniform canonical ensemble /Co 
is topologically always a circle, it evolves into the closed curve JC n after the n-th full 
period, with the preserved, constant, area enclosed by K n . This curve is just rotating 
and oscillating with n in case |Tr$i| = |5| < 2. If \S\ > 2, the action by $ n = $" is 
exponentially stretching in the direction of the eigenvector with eigenvalue ex > 1, and 
exponentially contracting in the direction of the other eigenvector with the eigenvalue 
C2 = 1/ei < 1. Therefore, the energy of the individual initial condition described by 
the formula (0) will be exponentially increasing for any initial condition (qo,Po) (vector 
in the phase plane (q,p)), except for the case when (qo,po) is exactly in the direction of 
the second eigenvector (stable manifold) corresponding to e 2 = 1/ei < 1. 

6. General formula for the energy evolution 

In this section we wish to consider an exact expression for the evolution of the 
energy distribution by studying a decomposition of one adiabatic process into several 
consecutive adiabatic processes. 

Let us first mention that the energy distribution P{E\) evolved from the original 
delta-like distribution 5(E — E ) is a kind of a Green function for the energy evolution. 
Let us denote it by G{E\] E ). If we have a distribution of initial energies w(E ), such 
that at each energy Eq the distribution on the energy contour is a uniform canonical 
distribution, then the final energy distribution is 



Thus by knowing G, which we call G-function, we can calculate the final energies of any 
family w(Eq) of uniform canonical ensembles of initial conditions. 

If the adiabatic process is ideal adiabatic, then the G-function is a delta function, 



For ensembles of other types, which are not uniform canonical, we must go back to 
our fundamental equation (J7J) and perform the averaging using the distribution in space 



Now suppose that the interval of length T is divided into an arbitrary number of 
finite subintervals (tj,tj + i), where t is the beginning of the process (interval) and t n is 
the end of the process, and j = 0, 1, . . . , n — 1. The behaviour of u(t) inside each j-th 
subinterval is assumed to be entirely arbitrary, but the process must be such that at 




(106) 



G{E 1 ]E ) = 6(E 1 -u 1 E /u ). 



(107) 



(Eo,<j>). 
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each integration step tj the distribution is uniform canonical. This condition is certainly 
satisfied if the process is ideal adiabatic, in general not. 

It is then obvious that the energy G-function G(E] E ) for the complete process 
divided into n subintervals is given by the multiple integral 



G(E;E )= (108) 
J ...J G n {E)x n -i)G n -i{x n -x]x n -^) . . .Gi{xi] Eo)dx n -i . . .dx2dxi. 



n-l 



Indeed, using (J1U7)) we can immediately verify this equation. All moments of the final 
distribution can be easily calculated as they are all fully determined by the first moment 
alone, so what we need is just the first moment. From our theory in section 2 we know 
that in full generality the first moment of any G(E; Eq) is a linear function of the initial 
value E , namely 

E = J EG(E; E )dE = gE , (109) 

where the constant g — (a + j3)/2 is a constant independent of E and is determined by 
the nature of u(t) inside the relevant interval of evolution. We shall call g the ^-factor of 
G. For an ideal adiabatic process we know from (|T2*|) that gj = uij/uj-x- From equation 
()108j) we immediately find the factorization property 

9 = 9n9n-i ■ ■ ■ 929u E = gE Q = g n . . .g 2 g 1 E Q . (110) 

Obviously, for an ideal adiabatic process where each gj = Uj/ujj-i, the above equation 
is certainly satisfied. 

It is possible also to show the converse [21] : If t ne composition formula ()108|) is 
true for any intermediate points of integration tj and Xj, then the process must be ideal 
adiabatic, implying that 

G{E f Ej-i) = 5(E 3 - UjEj.x/uj-x) (HI) 

applies for all j, and gj = Uj/ujj-x- This can be shown by splitting the time 
interval (to,t n ) into infinitesimal subintervals and using a piecewise constant function 
to approximate u(t), and then using gj = ^(^j/^j-i + -0 f rom equation for all j, 
finally evaluating g by equation (jllOj) . and finding g = uj n /ujQ, which implies that the 
process is ideal adiabatic at all times of the time interval, because /i 2 = 0. 

The composition formula ()108j) (factorization property of the G-function) will apply 
also in nonlinear systems, but the relationship between E\ and Eq is then no longer 
linear. Therefore using the composition formula ()108|) for infinitesimal intervals, and 
approximating u)(t) by piecewise constant or piecewise linear functions etc. might be 
of extreme importance to find new global powerful approximations for G-functions and 
their moments. 

The theory for nonlinear systems is left open for the future work. 
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7. Discussion and conclusions 

Our aim in this work is to study the time evolution of the energy in a general (no 
restriction upon u(t)) time- dependent ID harmonic oscillator in a rigorous way, and 
then also to calculate the statistical properties of the final energy distribution P{E\) 
after some time of lentgh T, if at the beginning we have a uniform canonical distribution 
of initial conditions at constant energy E . We are able to calculate rigorously all the 
moments of P(E\). Odd moments are exactly zero, the even moments are powers of the 
second moment and the second moment /i 2 is a function of the first moment. Therefore 
everything is determined by the first moment E\ and the variance /z 2 , which is a function 
of Ei. Thus the energy distribution function P{Ei) is universal in the sense that it does 
not depend on any other properties of u(t). (Mathematically speaking, the distribution 
function P{E\) is a two parameter family parametrized by E\ and /i 2 , but in physics 
these two parameters are connected through the formula (fTTj) . Of course, the shape of 
the distribution function P{E\) does not change if E\ is translationally shifted, and thus 
for the mathematical purposes we can set E\ = 0.) In this analysis we clearly see when 
the adiabatic invariant I(t) = E(t)/u)(t) is conserved or not. In the (ideal) adiabatic 
limit T —>■ oo it is conserved, the variance /i 2 is zero and E\ = E\ = uj\Eq/ujq. If it is 
not conserved exactly, when T is finite, we find /z 2 > 0, and it can be calculated using a 
WKB method analytically in a closed form, which is a major achievement of this work. 
We have also studied three specific solvable models and have demonstrated the power of 
the WKB expansion, where already the leading WKB term usually very well describes 
the asymptotic behaviour of /z 2 when e = 1/T goes to zero. We also show what happens 
if u)(t) is smooth and of class C m , having m continuous derivatives, calculating and 
showing that /j 2 oscillates as e goes to zero, but in the mean vanishes as oc e 2 ( m+1 ). If 
u(t) is analytic, thus it also is of class C°°, it is known from the literature that /x 2 must 
decay exponentially cx exp(const/e). If co{t) is periodic, E\ can grow exponentially, and 
so does the variance /i 2 , in which case l(t) = E{t)/ui{t) is not conserved, but we can 
describe the system. We have introduced the so-called G-function, which is a kind of 
a Green function for the evolution of the energy and derived a composition formula for 
it when the interval of evolution is decomposed into a finite number of subintervals, 
for which the corresponding G 3 -function is known for all subintervals j and is uniform 
canonical there. This formula applies also to nonlinear systems and might be a good 
starting basis to describe them. The theory for nonlinear systems remains open and is 
a subject of the current research 
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